Probability 2 (Data science Master, University of Lille) / CMF, (Centrale Lille, G3 SDIA)

---

Lab 1 - Discrete time homogeneous Markov chains

---

Guidelines (read carefully before starting)

Objectives: numerically simulate basic Markov chains (discrete time and discrete state space).

Setup: after retrieving the resources for the lab on moodle:

• place the .zip archive in a local folder (Computer -> Documents/Python/);
• unzip the archive .zip;
• rename the folder with the convention lab1_Name1_Name2
• duplicate the notebook file and rename it lab1_Name1_Name2.ipynb;
• [optional, possibly needed if working from Centrale's machines]
  • create a lab1 conda environment from the provided requirement.txt file

  conda create --name=lab1 --file=requirement.txt
  conda activate lab1
  # do not forget to deactivate the environment if needed
  # you can remove the environment once you are done
  conda env remove --name=lab1
  

  • launch jupyter notebook (the python environment will be the one from the conda environment lab1)
• at the end of the session, do not forget to transfer your work to your own network space if you are working on a machine from the school (do not leave your work on the C: drive).

Assessment → grade /20 (possibly converted later on to a grade ranging from F to A (A+))

This lab session will be evaluated, based on your answer to the exercises reported in a Jupyter notebook (e.g., this one) and any additional .py file produced. In particular:

• make sure the notebook you produce is properly annotated: each exercise should be introduced by a short sentence summarizing its context. Concisely answer each question from the guideline.
• relate the numerical results to the theory covered during the lecture whenever appropriate;
• codes without any explanations (in a text cell) will not be considered as a complete answer, and will only be attributed a fraction opf the grade attributed to the question.
• any code produced should be commented whenever appropriate;
• include appropriate axis labels and a relevant title to each figure reported in the notebook;
• document any function you introduce (using docstrings), and comment your code whenever appropriate (see, e.g., PEP 8 style recommendations).
  • use a reference docstring style, e.g., the google docstring style.
• give a clear structure to your notebook: longer theoretical explanations should be reported in markdown cells. Include LaTeX equations to improve the clarity of the document.

Additional evaluation criteria:

• Code clarity / conciseness / documentation
• Comment quality / overall presentation

Contents

• Exercise 1: Ehrenfest model
• Exercise 2: Simulation of a discrete time homogeneous Markov chain

%load_ext autoreload
%autoreload 2

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

from src.lab_1_utils import exercise_1_workflow, exercise_2_workflow

import plotly 

plotly.offline.init_notebook_mode()

---

Exercise 1: Ehrenfest model (↑)

Consider a system of $K = 30$ particles (labeled from 1 to $K$) evolving in a closed box. The box is divided into two compartments in contact with each other, respectively identified by an index, $0$ and $1$. A hole at the interface between the two compartments allows the particles to move from one compartment to the other.

Particle motion is modeled as follows: at each discrete time intant $n$, one particle is chosen uniformly at random and moved from its current compartment to the other. Let $X(n)$ denote the number of particles in compartment $0$ at time $n$.

1. Briefly justify that $\bigl(X(n) \bigr)_{n \in \mathbb{N}}$ is a Markov chain.

En décrivant le diagramme d'état de ce système on voit qu'il s'agit d'une chaîne de Markov :

Source : Computational Physics, Y. D. Chong, Nanyang Technological University

Dans notre cas, q=1 et N=30.

2. Is the chain irreducible? (Positive) recurrent? Aperiodic or periodic? Prove each statement from your answer (recall the main steps covered during the exercise session).

• Il n'y a qu'une seule classe de communication car il y a toujours un chemin de probabilité non nulle entre chaque état. Ainsi la chaîne de markov est irréductible.

• La chaîne de Markov est périodique de période 2 : En effet, si nous décomposons $S = \{0, ..., 30\}$ en deux sous parties : Les entiers pairs et les entiers impairs, nous passons forcément d'une sous partie à l'autre à chaque étape. C'est la plus petite décomposition possible, donc la chaîne est de période 2.

• Enfin, la chaîne est positive récurrente : Elle ne possède qu'une classe d'équivalence fermée.

3. Recall the structure of the transition matrix, and encode it in Python (without any for loop).

Hint: use the function numpy.diag.

On a, $ \forall i,j \in [|1, n|] $, $$ P(X_{n+1}=j | X_{n}=i) = \begin{cases} \frac{K-i}{K} &\text{si } j = i+1 \\ \frac{i}{K} &\text{si } j = i-1 \\ 0 &\text{sinon} \end{cases} $$

Initialisation du système et de la matrice de transition P :

K = 30
workflow_1 =  exercise_1_workflow(K)

workflow_1.initialize_P()

4. Numerically verify that the binomial distribution $\mathcal{B} (K, 1/2)$ is invariant for the chain. This invariant distribution will be denoted by $\pi$ in the rest of the exercise.

Hint: you can use scipy.stats.binom.

from scipy.stats import binom
binom_stat = binom(K, 0.5)

if workflow_1.check_invariant_distribution(binom_stat):
    print("Cette distribution est bien invariante pour notre système.")
else:
    print("Cette distribution n'est pas invariante pour notre système.")

Cette distribution est bien invariante pour notre système.

5. Implement a Python function ehrenfest to simulate a trajectory of the chain for a system of $K$ particles, for initial distribution $\mu$ (for instance a Dirac centered in $0$, meaning that the compartment $0$ is empty at $n = 0$). The maximum number of time steps will be controlled by an input parameter $n_{\max}$. Random number generation will be controlled by a random number generator passed as an input to the function.

For an efficient implementation, do not use vector-matrix product with the transition matrix implemented in 3: rely on the description of the system instead.

6. Simulate a trajectory of the chain starting in state $0$ for $n_{\max} = 5000$. Display the result in function of the time index $n$. Briefly describe the curve you obtained.

workflow_1.ehrenfest(5000)

On observe que le nombre de particules orbite autour de 15 passant d'abord au dessus de 15 puis en dessous. Cette simulation est cohérente avec la probabilité invariante trouvée. Pour en être sûre, nous pouvons regarder l'histogramme des états et le comparer a celui de la fonction binomiale :

workflow_1.plot_two_histograms()

Les deux histogrammes coïncident, nous obtenons bien une simulation de probabilité invariante $\pi$ = $\mathcal{B} (K, 1/2)$.

8. a) Modify the function defined in 1. so that it returns the return time to state 0, defined as $T_{0,0} = \inf \bigl\{ n > 0, X(n) = 0 \mid X(0) = 0 \bigr\}$.

T = workflow_1.average_return_time(n_average=1)

Pour K = 30, impossible de trouver un retour a 0.

En 100 000 steps il n'y a pas eu de retour à 0. Nous pouvons faire l'hyopothèse que pour 30 particules le temps de retour en 0 est difficile à observer car il peut prendre une valeur très grande, difficile a calculer en un temps raisonnable.

8. b) [Optional] Run several chains (about 5, ideally in parallel) for $K = 10$, $n_{\max} = 5000$, and compare the empirical average of $T_{0,0}$ to $\pi(0)$. What do you observe?

Hint: a good tutorial showing how to run functions in parallel in Python is available here.

workflow_1_10 = exercise_1_workflow(10)

T_10 = workflow_1_10.average_return_time()

Le temps de retour moyen en 0 pour 10 particules est de : 417.2 steps.

8. c) Comment on the possibility of numerically observing the chain returning to its initial state as $K$ increases.

workflow_1.average_return_time_evolution(n_average = 100)

D'après le graphique, le temps de retour moyen vers 0 semble suivre une courbe exponentielle en K. Il va donc très vite devenir difficile d'observer un retour en 0 en un temps de calcul raisonnable. Cependant, les états étant récurrents, il y aura toujours un retour en 0.

---

Exercise 2: Simulation of a discrete time homogeneous Markov chain (↑)

Let $\bigl( X(n) \bigr)_{n \in \mathbb{N}}$ be a discrete time homogeneous Markov chain defined by the following initial distribution $\mu$ and transition matrix $P$

$$ \mu = [0, 1, 0], % \quad % P = \begin{pmatrix} 0.2 & 0.7 & 0.1 \\ 0.9 & 0 & 0.1 \\ 0.2 & 0.8 & 0 \\ \end{pmatrix}. $$

1. What can you say about the Markov chain $X$? (irreducibility, positive recurrence, periodicity, ...). Justify each of your claim, citing the relevant results from the lecture.

image.png

Voici le diagramme d'état de cette chaîne de markov :

• Cette chaîne de Markov est irréductible : Il n'y a qu'une seule classe de communication. Comme pour l'exercice précédent, il y a toujours un chemin possible entre tous les états.

• Cette chaîne est apériodique : Pour chaque état, la probabilité d'aller dans n'importe quel autre état est non nulle. Par conséquent, la chaîne est forcément apériodique.

• Cette chaîne est positive récurrente : Elle ne possède qu'une classe d'équivalence fermée.

2. Write a function simulate_dthmc simulating the trajectory of the Markov chain $\bigl( X(n) \bigr)_{n \in \mathbb{N}}$ for $n_{\max}$ time steps. The signature of the function should include the following elements:

• list of inputs: transition matrix $P$, initial distribution $\mu$, number of time steps $n_{\max}$, random number generator;
• output: array of lenght $n_{\max}$ representing the trajectory.

To this end, do not use matrix vector products, which would lead to an extremely inefficient algorithm in this case.

Initialisation du système :

workflow_2 = exercise_2_workflow()

3. Simulate a trajectory of the chain for $n_{\max} = 2000$ starting from $X(0) = 1$. Plot the histogram of the states visited by the chain.

workflow_2.simulate_dthmc(n_max=2000)

4. Determine numerically an invariant distribution $\boldsymbol{\pi}$ of the chain (e.g., based on an eigendecomposition numpy.linalg.eig). Is it unique? Compare it to the histogram obtained in 2 (superimpose graphs). What can you conclude?

workflow_2.compute_invariant_distribution()

La distribution invariante est unique.
La distribution invariante est : [0.49197861 0.4171123  0.09090909]

Les vecteurs propres liés à la valeur propre 1 étant linéairement dépendants il n'existe qu'une seule distribution invariante pour une chaîne de markov : Le vecteur propre associé à la valeur propre 1 de norme 1.

workflow_2.plot_two_histograms()

Les historgrammes coïncident, nous pouvons faire l'hypothèse que cette distribution invariante est la distribution limite de cette chaîne de Markov.

5. a) Compute $\mu_n \triangleq \mu P^n$, the probability distribution of $X(n)$. What is the limit of $\mu_n$ as $n$ goes to $+\infty$? Illustrate the result numerically.

workflow_2.compute_probability_distribuction(100)

workflow_2.check_invariant_distribution()

La distribution de X(n) tend vers la distribution invariante obtenue à la question 5 : [0.49197861 0.4171123  0.09090909]

5. b) Display on the same graph the curves $n \mapsto \mu_n(i)$ for $i = 1, \dotsc , 3$, and compare with $\pi$. Display on another graph the function $n \mapsto \Vert \mu_n - \pi \Vert_1$, where $\Vert \cdot \Vert_1$ is the $\ell_1$ norm. What does each of these curves illustrate?

workflow_2.plot_mu_evolution()

On remarque que $\mu_{n}(i)$ converge très rapidement vers $\pi$. Le régime transitoire dure un peut moins de 20 steps puis le régime est stationnaire avec $\mu = \pi$

workflow_2.plot_norm_difference()

La convergence vers $\pi$ est très rapide !

6. For each state $i \in \{1, 2, 3 \}$, simulate 100 trajectories starting from the state $i$ until the return time to $i$. For each state, compute the (empirical) average return time. Compare with its theoretical value.

workflow_2.theorical_return_time()

Moyenne du temps de retour en 1 : 2.03
Moyenne du temps de retour en 2 : 2.4
Moyenne du temps de retour en 3 : 11.0

workflow_2.empirical_average_return_time()

Moyenne du temps de retour en 1 : 2.07
Moyenne du temps de retour en 2 : 2.43
Moyenne du temps de retour en 3 : 9.7

Pour 100 trajectoires, l'estimation n'est pas parfaite mais nous obtenons des valeurs proche de l'estimation théorique.

workflow_2.empirical_average_return_time(n_trajectories=10000)

Moyenne du temps de retour en 1 : 2.03
Moyenne du temps de retour en 2 : 2.4
Moyenne du temps de retour en 3 : 11.15

Pour 10 000 trajectoires le résultat est quasi celui de l'estimation théorique, les simulations confirment bien les résultats théoriques calculés.